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Abstract. In this paper we consider the mathematical model of thermo- and photo¬ 
acoustic tomography for the recovery of the initial condition of a wave field from 
knowledge of its boundary values. Unlike the free-space setting, we consider the wave 
problem in a region enclosed by a surface where an impedance boundary condition is 
imposed. This condition models the presence of physical boundaries such as interfaces 
or acoustic mirrors which reflect some of the wave energy back into the enclosed domain. 
By recognizing that the inverse problem is equivalent to a statement of boundary 
observability, we use control operators to prove the unique and stable recovery of the 
initial wave profile from knowledge of boundary measurements. Since our proof is 
constructive, we explicitly derive a solvable equation for the unknown initial condition. 
This equation can be solved numerically using the conjugate gradient method. We 
also propose an alternative approach based on the stabilization of waves. This leads 
to an exponentially and uniformly convergent Neumann series reconstruction when the 
impedance coefficient is not identically zero. In both cases, if well-known geometrical 
conditions are satisfied, our approaches are naturally suited for variable wave speed 
and for measurements on a subset of the boundary. 
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Submitted to: Inverse Problems 

1. Introduction 

Photoacoustic Tomography (PAT) and Thermoacoustic Tomography (TAT) are medical 
imaging modalities that combine the high contrast properties of electromagnetic 
waves, with the high resolution properties of acoustic waves. When a short 
electromagnetic pulse (laser/microwave) is used to irradiate a biological tissue, the 
photoacousic/thermoacoustic effect results in emission of acoustic signals that can be 
measured outside the object by wide-band ultrasonic transducers [34, 32, 5, 18, 31, 33]. 
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The goal of PAT and TAT is to obtain an image of the electromagnetic absorption 
properties of the tissne. These properties are highly related with the molecnlar structnre 
of the tissne, and hence conld reveal pathological conditions [36, 38, 26, 22], The process 
of recovering the actnal image of the tissne from the external data involves two steps. 
The hrst step consists of recovering the internal absorbed radiation from knowledge of 
the aconstic signal measured in the exterior. In both PAT and TAT, the hrst step is the 
same and is modeled as the reconstruction of an initial condition of a wave held from 
data measured on the boundary of the region. [5, 3]. 

In most of the mathematical analysis of TAT and PAT, there is the assumption 
that the acoustic waves propagate in open space [1, 14, 27, 29, 13, 30, 2]. However, in 
practice this is not always a valid assumption [35]. In some cases, the acoustic waves 
interact with boundaries/interfaces such as the skull, detectors, or the interface with air 
[17, 24, 21, 25, 16]. In other cases, certain boundaries are purposely introduced in the 
experimental setup to improve performance. For example, the bioengineering group at 
University College London observed that better imaging resolution and an artihcially 
increased detection aperture can be attained by partially or fully enclosing the target 
region with acoustic mirrors [7, 8, 10]. Another example was proposed in [15], where a 
45° angle acoustic rehector is used to improve the limited detection view of linear arrays. 
Improved image quality on phantom and ex-vivo experiments were observed. Moreover, 
their method can be implemented by having the acoustic rehectors on the boundary 
of an immersion tank. In such cases, the detectors will form an half-enclosed space 
for the object, enhancing angle detection and potentially improving the image quality 
substantially. In the mathematical literature this problem was hrst analyzed in [19, 12], 
where reconstruction algorithms are obtained for an isotropic medium in rectangular 
cavities. In [19] a method based on the fast Fourier transform was presented to study 
the case of constant speed of a rectangular domain. Holman and Kunyansky [12] used an 
eigenfunction expansion method to obtain an approximate reconstruction that improves 
as the measurement time increases. 

In our work, we recognize that this inverse problem is equivalent to a statement of 
observability in the theory of control for partial diherential equations (PDFs). Under 
the so-called geometric control condition (see Assumption 3.1), uniqueness and stability 
follows from the work of Bardos, Lebeau and Rauch [4]. This control approach is 
naturally suited for variable anisotropic wave speed and for measurements on a subset 
of the boundary. Also, to model the presence of a physical boundary, such as detectors 
or acoustic mirrors, we consider the problem in a region enclosed by a surface where 
an impedance boundary condition is imposed. The impedance coefficient is allowed 
to vary over all non-negative values. Zero values of the impedance coefficient model 
acoustically hard surfaces, intermediate values model absorption, and relatively high 
values approximate acoustically soft boundaries (such as the interface with air). 
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We present two different reconstrnction schemes for this problem. The hrst one is 
based on the proof of nniqueness and stability. We explicitly derive a solvable equation 
for the unknown initial condition. This equation can be solved using the conjugate 
gradient iterative method. The second approach is based on a Neumann series expansion 
introduced by Stefanov and Uhlmann [28] for the case of free space propagation. We show 
that the Neumann series reconstruction is valid when the geometric control condition 
holds for the region of the boundary where the impedance coefficient is positive. It is 
worth mentioning that this work was developed simultaneously to the work of Stefanov 
and Yang [37] where they present an average time reversal approach to obtain a Neumann 
series reconstruction in the case when the impedance coefficient is identically zero. 

The structure of this paper is as follows. Section 2 contains the mathematical 
formulation of the forward acoustic problem where we incorporate the impedance 
boundary condition and the possibility of anisotropic variable media modeled in 
geometric terms. Anisotropy is commonly encountered in biological media such as bones 
and muscles, and the Riemannian geometric formulation is a mathematically convenient 
manner to model this type of media [9, 28]. In Section 3 we briefly review some tools of 
control theory which we employ in the remainder of the paper. These tools are naturally 
suited for variable wave speed and for measurements on a subset of the boundary. Section 
4 contains the main abstract results concerning uniqueness and stability of the inversion, 
whereas Sections 5 and 6 describe our proposed constructive algorithms to recover the 
initial acoustic prohle from the boundary measurement. Finally, we provide a brief 
discussion of these results and of future work in Section 7. 

2. Mathematical Formulation 

In this paper we study the thermoacoustic tomography problem in the presence of 
heterogeneous media and an enclosing surface. More precisely we assume anisotropy 
being modeled in geometric terms. For dimension n > 2, let hi C M” be a smooth, 
bounded and simply connected domain with boundary dfl. Let g be a Riemannian 
metric dehned on hi. The propagation of acoustic waves is governed by the following 
system. 


d^u — Au = 0 

in 

(0, r) X H 

(la) 

u = Uq and dfU = 0 

on 

{t = 0} X H 

(lb) 

duU + Xdtu = 0 

on 

(0, r) X dQ 

(Ic) 


where uq G (hi) is the unknown initial acoustic profile. Here A denotes a second-order 
elliptic differential operator of the form 

Au = c^AgU — qu, 


( 2 ) 
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where c : —)■ M is the wave speed, Ag is the Laplace-Beltrami operator associated with 
the Riemannian metric g and g : —)■ R is a potential. Also, A : dfl —)■ R denotes the 
impedance coefficient at the bonndary dQ. We assnme that A is a non-negative fnnction. 
Physically, A = 0 models an acoustically hard surface and A —)■ cxd approximates an 
acoustically soft boundary. In general we allow A to vary on the boundary dQ to model 
the presence or absence of detectors, or the inhomogeneous nature of an enclosing surface. 
The length of the observation window of time is given by r < cxo which is dehned in 
the next section. In (Ic), the symbol dy denotes the conormal derivative with respect 
to the metric g at the boundary dVt. All the coefficients in (la)-(lc) are assumed to be 
sufficiently smooth. 

In order to consider partial measurements, we divide the boundary as the disjoint 
union dVL = P U \ P) and we assume that Kuq = M|(o,r)xr constitutes our 
measurements. Then, the goal of the thermoacoustic tomography problem is to hnd 
the initial prohle uq from knowledge of Auq- This is a challenging problem with intricate 
dependencies on the nature of the manifold (ff,c“^g), the partial boundary P and the 
time interval (0,t). We recognize that, in mathematical terms, the thermoacoustic 
tomography problem coincides with a problem of boundary observability, which is one of 
the central concepts of control theory for PDEs [11, 4, 9]. Hence, the objective of this 
paper is to constructively employ the tools of control theory to solve the thermoacoustic 
tomography problem in an enclosed domain with heterogeneous sound speed. 

3. Background on Control Theory 

Our approach relies on exact boundary controllability for the wave equation. Hence, the 
purpose of this section is to review some background and dehne the notation concerning 
control theory. For details, we refer the reader to [11, 4, 9, 20]. The hrst point is to 
ensure that the manifold (0,c“^g), on which the wave equation is dehned, yields exact 
boundary controllability. 

Following Bardos, Lebeau and Rauch [4], we assume that our problem enjoys the 
geometric control condition (GCC) for the Riemannian manifold (fl, c“^g) with only a 
portion P of the boundary dfl being accessible for control and observation. We assume 
that P is a smooth and simply connected domain relative to dfl. In this paper, we 
assume that P satishes the following. 

Assumption 3.1 (GCC). The geodesic how of (fl,c“^g) reaches P after possible 
rehections on dVt \ P in hnite time r, and {A > 0} C P. In other words, there exists 
T < oo such that any geodesic ray, originating from any point in at f = 0, eventually 
reaches P in a non-diffiactive manner (after possible geometrical rehections on dfl \ P) 
before time t = t. We refer to r as the exact controllability time. 

Throughout the paper {A > 0} denotes the set {x G dfl : X{x) > 0}, similarly 
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{A = 0} denotes the set {x G dfl : X{x) = 0}. The condition {A > 0} C T is a natural 
exact controllability condition for the initial problem (3 a)-(3(f) with Robin boundary 
conditions. Concerning the inverse problem, this condition guarantees uniqueness, 
stability and reconstruction using a conjugate gradient method (see Sections 4-5). 
However for a Neumann series reconstruction we need to assume that {A > 0} = T 
(see Section 6). 

We work with the standard Sobolev spaces based on square-integrable functions. 
The associated inner-product extends naturally as the duality pairing between 
functionals and functions. We should interpret the Hilbert space H^{Q) with the inner- 
product appropriately weighted by c“^| det g|^/^ so that c^Ag is formally self-adjoint with 
respect to the duality pairing of if°(f2). Sometimes this is written as if°(r2; c“^dVol), 
where dVol(x) = | detgl^A^i^;^ however to simplify notation and since there is no risk of 
confusion we will only write it as ff°(f2). The same weight, c“^| detgl^'^^, is incorporated 
in the inner product for This implies that the differential operator A is formally 

self-adjoint. 

Now we consider the following auxiliary problem. Given G if°((0,r) x T), find 
the generalized solution ^ G (^^([0, r]; for f = 0, 1 of the following problem 



in 

(0, r) X fl, 

(3a) 

^ = 0 and dA = 0 

on 

O 

X 

p 

(36) 

dA + Xdti = C 

on 

(0,r) X F, 

(3c) 

dA = 0 

on 

(0,r) X {dn\T). 

{Sd) 


Since we assume that A > 0 in the boundary conditions (3c) then the above is a well- 
posed problem satisfying a stability estimate of the following form, 

ll^llc'=([o,r];r^i-fc(o)) + l|gradg^||c-fc([o_r];/^-fc(n)) < ||CI|rfO((o,r)xr)5 
for k = 0,1. For details see [11, 4] and references therein. Here and in the rest of the 
paper, the symbol ;< means inequality up to a positive constant. This gives rise to the 
following definition of a bounded operator. 

Definition 3.2 (Solntion Operator). Let the solution operator 

5 ; X T) ^ 

be given by the map ( i—)■ dt^\t=r where ^ is the solution of (3 a)-(3d). 

The goal of the control problem is to prove that the operator S is surjective. 
In other words, given arbitrary cj) G the goal is to find a boundary condition 

c G id°((0,r) X T) to drive the solution ^ of (3 a)- (3d) from vanishing Cauchy data at 
t = 0 to the desired Cauchy data (^, dA) = (0, 0) at time t = r. For our own reference 
hereafter, we state the well-posedness of this control problem as a theorem. A proof 
obtained under the geometric control condition is found in [4]. See also the extension in 
[6] to less regular domains and coefficients. 


Multiwave imaging in an enclosure with variable wave speed 


6 


Theorem 3.3 (Exact Controllability). Let the geometric control condition 3.1 hold. 
Then for any function (j) G there is a boundary control C G r) x F) so that 

the solution ^ of (3a)-(3d) satisfies 

= (0) 0) time t = T. 

Among all such boundary controls there exists Cmin which is uniguely determined by 
cf) as the minimum norm control and satisfies the following stability condition 

||Cmin||j/0((0,r)xr) < C'H 

for some positive constant C = F, A, r). 

Along with the above theorem, we readily obtain the following dehnition of the 
control operator that we will nse in the rest of the paper. 

Definition 3.4 (Control Operator). Let the control operator 

C : if°((0,r) X F), 

be given by the map cf i—)■ C min where Cmin is dehned in Theorem 3.3. 

This leads to the following lemma which states the basic properties for the control 
operator. The proof is a direct conseqnence of Theorem 3.3. 

Lemma 3.5. The control operator given by Definition 3.4 is bounded. Moreover, we 
have that cf = SCcf, and the following estimates, 

||0||ho(o) < ||C0||HO((o,r)xr) < ll^llrfo(n), (4) 

hold for all cf G 

4. Uniqueness and Stability using Control Theory 

In this Section we state the inverse problem for thermoaconstic tomography and we also 
introdnce and prove onr first main resnlt. We start by recalling the forward problem 
(la)-(lc) which is well-posed with a solntion u G (^^([O, r]; for A: = 0,1 for an 

initial condition uq G H^{VL). 

Definition 4.1 (Inverse Problem). FTnder the geometric control condition 3.1, find 
the nnknown initial condition uq of the forward problem (la)-(lc) from knowledge of 
the partial Dirichlet data Awo = 'u|(o,T)xr- 

Now we are ready to rednce the inverse problem to a solvable eqnation using 
duality arguments. First, let us introduce an auxiliary time-reversal operator U : 
ifO((0,r) xD)^ Lf°((0,r) x D) given by 

{Uv){t, x) = v{t — t,x). 


( 5 ) 
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The operator U is clearly unitary and self-adjoint. The same is true for this operator 
dehned as U : i^°((0,r) x T) —)■ //°((0,r) x T). Also notice that UdtV = —dtUv for 
functions that posses a weak derivative in time. 

In what follows, we will evaluate the duality pairing between the terms in equation 
(la)-(lc) against xjj = where ^ is the solution of (3a)-(3d) with ( = C(p and (p G 
Now, since u is a weak solution of (la)-(lc) then it satishes, 

{dtu,'tjj)n\lzl = ^o^r)xdn - (gradgU, gradi/^)(o,r)xn 

- (o^r)xn + {dtU,dt'tlj)(o,T)xn- 

Similarly, ip = UP, is also a weak solution to the wave equation (3a), then it satishes, 

{u,dt'tp)n\lZl = {.u-,dy'ip){o,T)xda - (gradgW, gradi/^)(o,r)xo 

- (ri,c"^g'0)(o,r)xn + 9tt/')(o,r)xn- 

Subtracting the above identities and using the Cauchy data u = uq and dtu = 0 at 
f = 0, and the vanishing Cauchy data ip = dtip = 0 at time t = t, and the fact that 
dtip\t=o = —(p) we obtain, 

(itO, 0)f2 (0,T)xdfl (p,T)xdQ, 

= {u,d^ip - Xdtip)(o,r)xan - {duU + \dtu,ip)(o,T)xdn- 

Now notice that d,,ip = dJJp = Udyp and \dtip = XdtUp = —UXdtp on (0,r) x <917. 
Therefore, from the boundary conditions (Ic) and (3c)-(3d), we arrive at 

{uo,(p)n = {Auo,UC){o,r)xr = {Auo,UC(p)^o,r)xr, for all cp G H°{n), 

where Auq = u|(o,T)xr is the measurement at the partial boundary T. 

As a result, the unknown initial condition uq is explicitly recovered as follows, 

uq = C*UAuo, (6) 

where C* : if°((0,r) x T) ^ is the adjoint of C : —)■ H^((0,t) x T) which 

is a bounded operator. This proves the hrst of our main results. 

Theorem 4.2 (Uniqueness and Stability). The inverse problem of thermoacoustic 
tomography, as stated in Definition 4-fi is uniquely solvable. Moreover, the following 
stability estimate 

||wo||_h'0(q) < C'||AMo||HO((o,T)xr), 

holds for a positive constant C = C{Vt, T, A, r). 

Remark 4.3. From the details in [4], we have that the control operator C is bounded 
on Sobolev spaces of any scale of regularity (positive and negative). Therefore, we easily 
obtain from (6) a stability estimate of the form ||Mo||jri(o) ^ ||AMo||rri((o,T)xr) which 
matches the level of regularity assumed in the forward problem (la)-(lc). 
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In order to take full advantage of the recovery equation (6), we now proceed to express 
the adjoint-control operator C* in terms of applicable operators. Hence, we remind the 
reader of the following facts from control theory. For a proof see [11, Ch. 6] and [4, 20]. 

Lemma 5.1. Let the mapping S* : —)■ if°((0,r) x F) be given by z ^ 'W^|(o,T)xr 

where w G (^^([O, rj; for /c = 0,1 is the generalized solution of the following 
time-reversed problem, 


d^w — Aw = 0 

in (0, t) X Q 

(7o) 

w = z and dtw = 0 

on {t = t} X D 

(76) 

d,yW — XdfW = 0 

on (0, r) X dQ. 

(7c) 

notation suggests, S* 

is the adjoint of the solution 

operator S from 


Definition 3.2. Moreover, if the geometrical control condition 3.1 holds, then we have: 


1. The operator (SS*) : —)■ is coercive. 

2. The control operator C is given by C = S*{SS*)~^. 

3. An eguation of the type {SS*)(j) = (3 can be solved using the conjugate gradient 
method. 


Our second main result follow directly from Lemma 5.1 and the explicit formula for 
Uq obtained in (6). 

Theorem 5.2 (Reconstruction Method 1). The adjoint of the control operator C 
satisfies C* = {SS*)~^S. Hence, the unigue solution Uq of the thermoacoustic tomography 
problem is given by, 

uo = {SS*)-^SUAuo, ( 8 ) 

where U is the time-reversal operator defined in (5) and (SS*) can be inverted using the 
conjugate gradient iterative method with convergence in 


6. Reconstruction by Neumann Series 

In this section we consider an alternative reconstruction using a Neumann series 
expansion in the case when F = {A > 0} is non empty and satishes the GCC 3.1. 
In other words, we show that if F satishes the GCC 3.1 with {A > 0} C F replaced by 
F = {A > 0}, then we can recover uq from a Neumann series applied to a backprojection 
of the data restricted to F. Our analysis follows mainly the stabilization of waves in the 
work by Bardos, Lebeau and Rauch [4, Sect. 5] and uses in a non-trivial way the energy 
decay guaranteed by {A > 0} = F. For the case A = 0 the construction presented here 
will not work, even in the case of full data, because a certain error operator K will fail 
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to be a contraction due to the fact that total energy in conserved. For this case, we refer 
the reader to [37] were an average time reversal approach is used to avert this difficulty 
and obtain a Neumann series expansion when A = 0. 

We use a similar backprojection operator as the one introduced by Stefanov and 
Uhlmann [28] for the free-space setting. See also [23]. However in our case, we can 
use the information of the Neumann derivative of the wave in {A > 0} to make the 
error operator K satisfy Neumann boundary condition in all of dVt. This avoids using 
a harmonic extension at f = r and allows us to obtain a Neumann series reconstruction 
with information only on the partial boundary F. 

We introduce some notation for energy spaces with Robin boundary condition and 
describe some properties of the evolution operator for the Cauchy problem. 

6.1. Energy Spaces and Evolution Operators 
Consider the initial Cauchy value problem 


dtu 

— Au = 0 

in 

(0, r) X D 

(9 a) 

u = 

Mo and dtU = ui 

on 

{t = 0} X D 

{9b) 

dyU 

-|- XdtU = 0 

on 

(0, r) X dOt 

(9c) 


for A a non-negative smooth function. For this subsection we allow A = 0. 

Dehne the energy of u at time t by 

E{t, u) = \- [ IgradgUp -f c~‘^q\u\^ + |^dVol. (10) 

^ Jn 

Denote by CD^(r2) the set of Cauchy data related to Robin boundary condition 
(uoWi) foi' which there exist a unique smooth solutions of the initial problem (9a)-(9c) 
in (0, r) X D. 

Denote by CD)j(D) the closure of CD^(D) in i7^(D)xL^(D). In CD)j(D) we consider 
the energy norm 

ll(wo,Mi)||| = ^ / IgradgMo^ + + c^Vi^dVol. (11) 

^ Jn 

When q is not identically zero, the energy norm is equivalent to the H^{Q) x L^(D) 
norm in CD)j(r2). However, when q is identically zero, (11) is not a norm in CD)^(r2) 
since constants are non-trivial solutions of Au = 0, d^u + Xdtu = 0, and hence 
II (const., 0)||£; = 0. In such case, we consider the Cauchy data set CDJ^°(D), given 
by the {uo,Ui) G CD^(r2) such that 
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where dS is the induced Riemannian metric on with the induced orientation. It 
follows easily that (11) is also a norm on CD)^°(f2). To justify this condition (12), notice 


Ut{t,x)dx+ / A(a;)'u(f, a:)dS'= independent of time, 

Jan 


(13) 


hence by replacing u in (9a)-(9c), with 

In dx + \u dS 


u 


we reduce the non-uniqueness problem in (9 a)-(9c) to the study of solution for which 


/ Ut{t,x)dx+ / A(x)n(f, a;)dS = 0 
Jn Jan 

In summary, we consider the problem (9a)-(9c) with Cauchy data set CD)j(r2) when q 
is not identically zero and CDj^°(f2) otherwise; to simplify notation we denote this set 
by CD/j(n). Recall that in CD/^(f2) we always consider the norm || • H® given by (11), 
whether or not q is equals to zero. 

For the problem of thermoacoustic and photoacoustic tomography we deal only 
with the case Ui = 0 in (9a)-(9c). For this reason we denote by the projection of 

CD/j(f2) onto its hrst component. The norm induced in Hji{Q) is the Dirichlet norm 

ll“o||^« = ^ / IgradgWoP + c"^g|MoPdVol. (14) 

^ Jn 

Clearly is a subspace of CDjj(f2), moreover CD/j(r2) = H(B. The later 

follows from the fact that the energy of u at time t is given by E(t, u) = 

We now dehne the evolution operators related to Robin and Neumann boundary 
conditions. Let f > 0. First consider the case when {A > 0} is not empty. In this case, 
we denote by S^it) : CD 7 ^(f 2 ) —)■ CD 7 ^(A 2 ) the evolution operator of the wave equation 
related to Robin boundary conditions. The operator takes Cauchy data at time 

0, to Cauchy data at time t. Under the geometric control condition 3.1 with exact 
controllability time r, this evolution operator S'ij(f) dehnes a contraction semigroup in 
(CDij(U), II • ||e), see Theorems 5.5 in [4] for the case q not zero and Theorem 5.6 in [4] 
for the case q equals to zero. Notice that since 

dtE{t,u) = — f X{x)\ut\‘^dS, 

Jan 

Robin boundary conditions are energy decreasing only if {A > 0} is non empty. 
For the case when A = 0, denote by Sisfit) : CDAr(fl) —)■ CDAr(f2) the evolution 
operator of the wave equation related to Neumann boundary conditions. Since Neumann 
boundary conditions are energy preserving then ||5'Ar(t)|| = 1 as an isomorphism in 

{CI)^{Q),\\-\\e). 

We now present the Neumann series reconstruction for the case of full data. 
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Let M be a solution of the initial value problem (9a)-(9c). Let r > 0 be the exact 
controllability time, we dehne the “back projection operator” A : iL^((0,r) x dVL) —)■ 


H^iM) that acts on the data Kuq = n|(o,T-)x( 9 o by 
(74Ano)(ic) = n(0, x) 

where v solves 

d^v — Av = 0 in (0,r) x (15a) 

n = 0 and dtv = 0 on {t = r} x {15b) 

dyV = —XdtAuo on (0,r) x L (15c) 

dyV = 0 on (0, r) x dQ \ T (15d) 


We think of AAuo as the hrst approximation on the reconstruction of uq. Notice that 
we only need the data Auq on L. The “error operator” K is given by 

K := Id- AA 

Notice that if we let tc = n — n, then {Kuo){x) = ta(0,a;); moreover w solves 


d^w — Aw = 0 

in (0, t) X Q 

(16a) 

w = u{t, •) and dfW = dtu{T, •) 

on {t = t} X Q 

{16b) 

dyW = 0 

on (0, r) X dH 

(16 c) 


Then we can write K : Hji{Q) —)■ Hji{Q) as a composition of operators 

K = niSN{-r)SR{T)nl 

where 7r*(Mo) = (wo,0) and 7ii{uo,ui) = uq. 

The important observation is that, if T = {A > 0} satishes the GCC 3.1 with exact 
controllability time r, K is the composition of unitary operators and Sr{t) which is 
a contraction, see Theorems 5.5 and 5.6 in [4]. Hence K itself is a contraction, i.e., 
WKWhu^Hh < 1- This implies in particular that we can reconstruct Uq from the data 
Auolr via a Neumann series expansion. We state this observation as a theorem. 

Theorem 6.1 (Reconstruction Method 2). Let T = {A > 0} &e non empty and 
satisfying the GCC 3.1 with exact controllability time r. Then AA = Id — K, where K 
is a contraction, \\K\\Hj^^Hn < 1- In particular, Id — K is invertible on Hr{Q) and the 
inverse thermoacoustic problem has an explicit Neumann series reconstruction given by 

OO 

uo = Y^ K^AAuo 

n=0 


( 17 ) 
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7. Conclusion 

We have proposed two reconstruction algorithms for the thermoacoustic tomography 
problem in an enclosure with observation on a portion V of the enclosing boundary 
dVt. In both cases, the reflecting portion 9fl\r of the boundary allows the energy of the 
unknown initial condition to eventually reach V where measurements take place. Our 
work provides precise conditions under which such a scenario yields the solvability of this 
inverse problem, and explains from a mathematical viewpoint how the use of acoustic 
mirrors may render a good reconstruction for limited-view detection. 

The hrst algorithm (Theorem 5.2), based on the conjugate gradient method, is 
naturally suited for full enclosure where the acoustic energy may not be allowed to 
escape from the domain of interest. The second algorithm (Theorem 6.1), based on a 
Neumann series expansion, is computationally more attractive and naturally suited for 
partial enclosure where the energy is allowed to leave through the observation part of 
the boundary. The authors are in the process of developing numerical implementations 
to compare the proposed reconstruction methods with each other and with experimental 
data. As soon as meaningful results are obtained from these efforts, they will be reported 
in a forthcoming publication. 
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